Cours n°01 : Introduction à l’inférence bayésienne Cours n°02 : Modèle Beta-Binomial Cours n°03 : Introduction à brms, modèle de régression linéaire Cours n°04 : Modèle de régression linéaire (suite) Cours n°05 : Markov Chain Monte Carlo Cours n°06 : Modèle linéaire généralisé Cours n°07 : Comparaison de modèles Cours n°08 : Modèles multi-niveaux Cours n°09 : Modèles multi-niveaux généralisés Cours n°10 : Data Hackathon
Contraintes : Certaines valeurs soient fortement probables (autour de la moyenne \(\mu\)). Plus on s’éloigne, moins les valeurs sont probables (en suivant une décroissance exponentielle).
D’où vient la loi normale ?
\[
y = \exp \big[-x^{2} \big]
\]
On étend notre fonction aux valeurs négatives.
D’où vient la loi normale ?
\[
y = \exp \big[-x^{2} \big]
\]
Les points d’inflection nous donnent une bonne indication de là où la plupart des valeurs se trouvent (i.e., entre les points d’inflection). Les pics de la dérivée nous montrent les points d’inflection.
D’où vient la loi normale ?
\[
y = \exp \bigg [- \frac{1}{2} x^{2} \bigg]
\]
Ensuite on standardise la distribution de manière à ce que les deux points d’inflection se trouvent à \(x = -1\) et \(x = 1\).
Mais… cette distribution n’intègre pas à 1. On divise donc par une constante de normalisation (la partie gauche), afin d’obtenir une distribution de probabilité.
Modèle gaussien
Nous allons construire un modèle de régression, mais avant d’ajouter un prédicteur, essayons de modéliser la distribution des tailles.
On cherche à savoir quel est le modèle (la distribution) qui décrit le mieux la répartition des tailles. On va donc explorer toutes les combinaisons possibles de \(\mu\) et \(\sigma\) et les classer par leurs probabilités respectives.
Notre but, une fois encore, est de décrire la distribution postérieure, qui sera donc d’une certaine manière une distribution de distributions.
Modèle gaussien
On définit \(p(\mu, \sigma)\), la distribution a priori conjointe de tous les paramètres du modèle. On peut spécifier ces priors indépendamment pour chaque paramètre, sachant que \(p(\mu, \sigma) = p(\mu) p(\sigma)\).
On définit \(p(\mu, \sigma)\), la distribution a priori conjointe de tous les paramètres du modèle. On peut spécifier ces priors indépendamment pour chaque paramètre, sachant que \(p(\mu, \sigma) = p(\mu) p(\sigma)\).
Il s’agit de la même formule vue lors des cours 1 et 2, mais cette fois en considérant qu’il existe plusieurs observations de taille (\(h_{i}\)), et deux paramètres à estimer : \(\mu\) et \(\sigma\).
Pour calculer la vraisemblance marginale (en vert), il faut donc intégrer sur deux paramètres : \(\mu\) et \(\sigma\). On réalise ici encore que la probabilité a posteriori est proportionnelle au produit de la vraisemblance et du prior.
Distribution postérieure - Grid approximation
# définit une grille de valeurs possibles pour mu et sigmamu.list <-seq(from =140, to =160, length.out =200)sigma.list <-seq(from =4, to =9, length.out =200)# étend la grille en deux dimensions (chaque combinaison de mu et sigma)post <-expand.grid(mu = mu.list, sigma = sigma.list)# calcul de la log-vraisemblance (pour chaque couple de mu et sigma)post$LL <-sapply(1:nrow(post),function(i) sum(dnorm( d2$height,mean = post$mu[i],sd = post$sigma[i],log =TRUE) ) )# calcul de la probabilité a posteriori (non normalisée)post$prod <- post$LL +dnorm(post$mu, 178, 20, log =TRUE) +dunif(post$sigma, 0, 50, log =TRUE)# on "annule" le log en avec exp() et on standardise par la valeur maximalepost$prob <-exp(post$prod -max(post$prod) )
Distribution postérieure - Grid approximation
# select random 20 rows of the dataframe post %>%slice_sample(n =20, replace =FALSE)
Under the hood : Stan est un langage de programmation probabiliste écrit en C++, et qui implémente plusieurs algorithmes de MCMC : HMC, NUTS, L-BFGS…
data { int<lower=0> J; // number of schools real y[J]; // estimated treatment effects real<lower=0> sigma[J]; // s.e. of effect estimates }parameters { real mu; real<lower=0> tau; real eta[J]; }transformed parameters { real theta[J];for (j in1:J) theta[j] = mu + tau * eta[j]; }model { target +=normal_lpdf(eta |0, 1); target +=normal_lpdf(y | theta, sigma); }
Bayesian regression models using Stan
Le package brms(Bürkner, 2017) permet de fitter des modèles multi-niveaux (ou pas) linéaires (ou pas) bayésiens en Stan mais en utilisant la syntaxe de lme4.
se spécifie avec brms (comme avec lme4) de la manière suivante :
model <-brm(y ~ x + (1| subject) + (1| item), data = d, family =gaussian() )
Rappels de syntaxe
Le package brms utilise la même syntaxe que les fonctions de base R (comme lm) ou que le package lme4.
Reaction ~ Days + (1+ Days | Subject)
La partie gauche représente notre variable dépendante (ou outcome, i.e., ce qu’on essaye de prédire). Le package brms permet également de fitter des modèles multivariés (plusieurs outcomes) en les combinant avec mvbind().
mvbind(Reaction, Memory) ~ Days + (1+ Days | Subject)
La partie droite permet de définir les prédicteurs. L’intercept est généralement implicite, de sorte que les deux écritures ci-dessous sont équivalentes.
mvbind(Reaction, Memory) ~ Days + (1+ Days | Subject)mvbind(Reaction, Memory) ~1+ Days + (1+ Days | Subject)
Rappels de syntaxe
Si l’on veut fitter un modèle sans intercept (why not), il faut le spécifier explicitement comme ci-dessous.
mvbind(Reaction, Memory) ~0+ Days + (1+ Days | Subject)
Par défaut brms postule une vraisemblance gaussienne. Ce postulat peut être changé facilement en spécifiant la vraisemblance souhaitée via l’argument family.
brm(Reaction ~1+ Days + (1+ Days | Subject), family =lognormal() )
Lisez la documentation (c’est très enthousiasmant à lire) accessible via ?brm.
Quelques fonctions utiles
# générer le code du modèle en Stanmake_stancode(formula, ...)stancode(fit)# définir les priorsget_prior(formula, ...)set_prior(prior, ...)# récupérer les prédiction du modèlefitted(fit, ...)predict(fit, ...)conditional_effects(fit, ...)# posterior predictive checkingpp_check(fit, ...)# comparaison de modèlesloo(fit1, fit2, ...)bayes_factor(fit1, fit2, ...)model_weights(fit1, fit2, ...)# test d'hypothèsehypothesis(fit, hypothesis, ...)
Un premier exemple
library(brms)mod1 <-brm(height ~1, data = d2)
posterior_summary(mod1, pars =c("^b_", "sigma"), probs =c(0.025, 0.975) )
Ces données représentent les distributions marginales de chaque paramètre. En d’autres termes, la probabilité de chaque valeur de \(\mu\), après avoir moyenné sur toutes les valeurs possible de \(\sigma\), est décrite par une distribution gaussienne avec une moyenne de \(154.59\) et un écart type de \(0.41\). L’intervalle de crédibilité (\(\neq\) intervalle de confiance) nous indique les 95% valeurs de \(\mu\) ou \(\sigma\) les plus probables (sachant les données et les priors).
En utilisant notre prior
Par défaut brms utilise un prior très peu informatif centré sur la valeur moyenne de la variable mesurée. On peut donc affiner l’estimation réalisée par ce modèle en utilisant nos connaissances sur la distribution habituelle des tailles chez les humains.
La fonction get_prior() permet de visualiser une liste des priors par défaut ainsi que de tous les priors qu’on peut spécifier, sachant une certaine formule (i.e., une manière d’écrire notre modèle) et un jeu de données.
priors <-c(prior(normal(178, 20), class = Intercept),prior(exponential(0.01), class = sigma) )mod2 <-brm( height ~1,prior = priors,family =gaussian(),data = d2 )
En utilisant notre prior
summary(mod2)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: height ~ 1
Data: d2 (Number of observations: 352)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 154.61 0.42 153.79 155.44 1.00 3367 2534
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 7.76 0.29 7.23 8.37 1.00 3235 2603
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
En utilisant un prior plus informatif
priors <-c(prior(normal(178, 0.1), class = Intercept),prior(exponential(0.01), class = sigma) )mod3 <-brm( height ~1,prior = priors,family =gaussian(),data = d2 )
En utilisant un prior plus informatif
summary(mod3)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: height ~ 1
Data: d2 (Number of observations: 352)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 177.87 0.10 177.68 178.06 1.00 2529 2389
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 24.60 0.95 22.84 26.52 1.00 3041 2444
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
On remarque que la valeur estimée pour \(\mu\) n’a presque pas “bougée” du prior…mais on remarque également que la valeur estimée pour \(\sigma\) a largement augmentée. Nous avons dit au modèle que nous étions assez certain de notre valeur de \(\mu\), le modèle s’est ensuite “adapté”, ce qui explique la valeur de \(\sigma\)…
Précision du prior (heuristique)
Le prior peut généralement être considéré comme un posterior obtenu sur des données antérieures.
On sait que le \(\sigma\) d’un posterior gaussien nous est donné par la formule :
\[\sigma_{post} = 1 / \sqrt{n}\]
Qui implique une quantité de données\(n = 1 / \sigma^2_{post}\). Notre prior avait un \(\sigma = 0.1\), ce qui donne \(n = 1 / 0.1^2 = 100\).
On peut donc considérer que le prior \(\mu \sim \mathrm{Normal}(178, 0.1)\) est équivalent au cas dans lequel nous aurions observé \(100\) tailles de moyenne \(178\).
Visualiser les échantillons de la distribution postérieure
post <-as_draws_df(x = mod2) %>%mutate(density =get_density(b_Intercept, sigma, n =1e2) )ggplot(post, aes(x = b_Intercept, y = sigma, color = density) ) +geom_point(size =2, alpha =0.5, show.legend =FALSE) +labs(x =expression(mu), y =expression(sigma) ) + viridis::scale_color_viridis()
Récupérer les échantillons de la distribution postérieure
On considère un modèle de régression linéaire avec un seul prédicteur, une pente, un intercept, et des résidus distribués selon une loi normale. La notation :
Les notations ci-dessus sont équivalentes, mais la dernière est plus flexible, et nous permettra par la suite de l’étendre plus simplement aux modèles multi-niveaux.
Dans ce modèle \(\mu\) n’est plus un paramètre à estimer (car \(\mu\) est déterminé par \(\alpha\) et \(\beta\)). À la place, nous allons estimer \(\alpha\) et \(\beta\).
Rappels : \(\alpha\) est l’intercept, c’est à dire la taille attendue, lorsque le poids est égal à \(0\). \(\beta\) est la pente, c’est à dire le changement de taille attendu quand le poids augmente d’une unité.
Régression linéaire à un prédicteur continu
priors <-c(prior(normal(178, 20), class = Intercept),prior(normal(0, 10), class = b),prior(exponential(0.01), class = sigma) )mod4 <-brm( height ~1+ weight,prior = priors,family =gaussian(),data = d2 )
Représenter l’incertitude sur \(\mu\) via fitted()
# on crée un vecteur de valeurs possibles pour "weight"weight.seq <-data.frame(weight =seq(from =25, to =70, by =1) )# on récupère les prédictions du modèle pour ces valeurs de poidsmu <-data.frame(fitted(mod4, newdata = weight.seq) ) %>%bind_cols(weight.seq)# on affiche les 10 premières lignes de muhead(mu, 10)
Pour rappel, voici notre modèle : \(h_{i} \sim \mathrm{Normal}(\alpha + \beta x_{i}, \sigma)\). Pour l’instant, on a seulement représenté les prédictions pour \(\mu\). Comment incorporer \(\sigma\) dans nos prédictions ?
# on crée un vecteur de valeurs possibles pour "weight"weight.seq <-data.frame(weight =seq(from =25, to =70, by =1) )# on récupère les prédictions du modèle pour ces valeurs de poidspred_height <-data.frame(predict(mod4, newdata = weight.seq) ) %>%bind_cols(weight.seq)# on affiche les 10 premières lignes de pred_heighthead(pred_height, 10)
Deux sources d’incertitude dans le modèle : incertitude concernant l’estimation de la valeur des paramètres mais également concernant le processus d’échantillonnage.
Incertitude épistémique : La distribution a posteriori ordonne toutes les combinaisons possibles des valeurs des paramètres selon leurs plausibilités relatives.
Incertitude aléatoire : La distribution des données simulées est elle, une distribution qui contient de l’incertitude liée à un processus d’échantillonnage (i.e., générer des données à partir d’une gaussienne).
d %>%# on utilise d au lieu de d2ggplot(aes(x = weight, y = height) ) +geom_point(colour ="white", fill ="black", pch =21, size =3, alpha =0.8)
Scores standardisés
d <- d %>%mutate(weight.s = (weight -mean(weight) ) /sd(weight) )d %>%ggplot(aes(x = weight.s, y = height) ) +geom_point(colour ="white", fill ="black", pch =21, size =3, alpha =0.8)
c(mean(d$weight.s), sd(d$weight.s) )
[1] -2.712698e-18 1.000000e+00
Scores standardisés
Pourquoi standardiser les prédicteurs ?
Interprétation. Permet de comparer les coefficients de plusieurs prédicteurs. Un changement d’un écart-type du prédicteur correspond à un changement d’un écart-type sur la réponse (si la réponse est aussi standardisée).
Fitting. Quand les prédicteurs contiennent de grandes valeurs, cela peut poser des problèmes de convergence (cf. Cours n°05)…
À vous de construire et fitter ce modèle en utilisant brms::brm().
Modèle de régression polynomiale
priors <-c(prior(normal(156, 100), class = Intercept),prior(normal(0, 10), class = b),prior(exponential(0.01), class = sigma) )mod6 <-brm(# NB: polynomials should be written with the I() function... height ~1+ weight.s +I(weight.s^2),prior = priors,family =gaussian(),data = d )
Modèle de régression polynomiale
summary(mod6)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: height ~ 1 + weight.s + I(weight.s^2)
Data: d (Number of observations: 544)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 146.66 0.37 145.93 147.39 1.00 3288 3142
weight.s 21.39 0.29 20.82 21.97 1.00 3027 2439
Iweight.sE2 -8.42 0.28 -8.98 -7.86 1.00 3578 2932
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 5.78 0.18 5.45 6.15 1.00 4048 2836
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Représenter les prédictions du modèle
# on crée un vecteur de valeurs possibles pour "weight"weight.seq <-data.frame(weight.s =seq(from =-2.5, to =2.5, length.out =50) )# on récupère les prédictions du modèle pour ces valeurs de poidsmu <-data.frame(fitted(mod6, newdata = weight.seq) ) %>%bind_cols(weight.seq)pred_height <-data.frame(predict(mod6, newdata = weight.seq) ) %>%bind_cols(weight.seq)# on affiche les 10 premières lignes de pred_heighthead(pred_height, 10)
Plusieurs méthodes pour calculer les tailles d’effet dans les modèles bayésiens. Gelman & Pardoe (2006) proposent une méthode pour calculer un \(R^{2}\) basé sur l’échantillon.
Marsman & Wagenmakers (2017) et Marsman et al. (2019) généralisent des méthodes existantes pour calculer un \(\rho^{2}\) pour les designs de type ANOVA (i.e., avec prédicteurs catégoriels), qui représente une estimation de la taille d’effet dans la population, et non basé sur l’échantillon.
Similar to most of the ES measures that have been proposed for the ANOVA model, the squared multiple correlation coefficient \(\rho^{2}\) […] is a so-called proportional reduction in error measure (PRE). In general, a PRE measure expresses the proportion of the variance in an outcome \(y\) that is attributed to the independent variables \(x\)(Marsman et al., 2019).
On a présenté un nouveau modèle à deux puis trois paramètres : le modèle gaussien, puis la régression linéaire gaussienne, permettant de mettre en relation deux variables continues.
Comme précédemment, le théorème de Bayes est utilisé pour mettre à jour nos connaissances a priori quant à la valeur des paramètres en une connaissance a posteriori, synthèse entre nos priors et l’information contenue dans les données.
La package brms permet de fitter toutes sortes de modèles avec une syntaxe similaire à celle utilisée par lm().
La fonction fitted() permet de récupérer les prédictions d’un modèle fitté avec brms (i.e., un modèle de classe brmsfit).
La fonction predict() permet de simuler des données à partir d’un modèle fitté avec brms.
Travaux pratiques - 1/2
Sélectionner toutes les lignes du jeu de données howell correspondant à des individus mineurs (age < 18). Cela devrait résulter en une dataframe de 192 lignes.
Fitter un modèle de régression linéaire en utilisant la fonction brms::brm(). Reporter et interpréter les estimations de ce modèle. Pour une augmentation de 10 unités de weight, quelle augmentation de taille (height) le modèle prédit-il ?
Faire un plot des données brutes avec le poids sur l’axe des abscisses et la taille sur l’axe des ordonnées. Surimposer la droite de régression du modèle et un intervalle de crédibilité à 89% pour la moyenne. Ajouter un intervalle de crédibilité à 89% pour les tailles prédites.
Que pensez-vous du fit du modèle ? Quelles conditions d’application du modèle seriez-vous prêt.e.s à changer, afin d’améliorer le modèle ?
Travaux pratiques - 2/2
Imaginons que vous ayez consulté une collègue experte en allométrie (i.e., les phénomènes de croissance différentielle d’organes) et que cette dernière vous explique que ça ne fait aucun sens de modéliser la relation entre le poids et la taille… alors qu’on sait que c’est le logarithme du poids qui est relié à la taille !
Modéliser alors la relation entre la taille (cm) et le log du poids (log-kg). Utiliser la dataframe howell en entier (les 544 lignes). Fitter le modèle suivant en utilisant brms::brm().
Où \(h_{i}\) est la taille de l’individu \(i\) et \(w_{i}\) le poids de l’individu \(i\). La fonction pour calculer le log en R est simplement log(). Est-ce que vous savez interpréter les résultats ? Indice: faire un plot des données brutes et surimposer les prédictions du modèle…
Proposition de solution
# on garde seulement les individus ayant moins de 18 ansd <-open_data(howell) %>%filter(age <18)priors <-c(prior(normal(150, 100), class = Intercept),prior(normal(0, 10), class = b),prior(exponential(0.01), class = sigma) )mod7 <-brm( height ~1+ weight,prior = priors,family =gaussian(),data = d )
Proposition de solution
summary(mod7, prob =0.89)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: height ~ 1 + weight
Data: d (Number of observations: 192)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-89% CI u-89% CI Rhat Bulk_ESS Tail_ESS
Intercept 58.23 1.40 56.00 60.46 1.00 4163 3249
weight 2.72 0.07 2.61 2.83 1.00 3866 3171
Family Specific Parameters:
Estimate Est.Error l-89% CI u-89% CI Rhat Bulk_ESS Tail_ESS
sigma 8.53 0.43 7.89 9.24 1.00 4146 2878
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Représenter les prédictions du modèle
# on crée un vecteur de valeurs possibles pour "weight"weight.seq <-data.frame(weight =seq(from =5, to =45, length.out =1e2) )# on récupère les prédictions du modèle pour ces valeurs de poidsmu <-data.frame(fitted(mod7, newdata = weight.seq, probs =c(0.055, 0.945) ) ) %>%bind_cols(weight.seq)pred_height <-data.frame(predict(mod7, newdata = weight.seq, probs =c(0.055, 0.945) ) ) %>%bind_cols(weight.seq)# on affiche les 6 premières lignes de pred_heighthead(pred_height)
# on considère maintenant tous les individusd <-open_data(howell)mod8 <-brm(# on prédit la taille par le logarithme du poids height ~1+log(weight),prior = priors,family =gaussian(),data = d )
Proposition de solution
summary(mod8, prob =0.89)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: height ~ 1 + log(weight)
Data: d (Number of observations: 544)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-89% CI u-89% CI Rhat Bulk_ESS Tail_ESS
Intercept -23.62 1.38 -25.84 -21.38 1.00 4071 3010
logweight 47.03 0.40 46.40 47.66 1.00 4120 2934
Family Specific Parameters:
Estimate Est.Error l-89% CI u-89% CI Rhat Bulk_ESS Tail_ESS
sigma 5.16 0.16 4.92 5.41 1.00 4314 3248
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Représenter les prédictions du modèle
# on crée un vecteur de valeurs possibles pour "weight"weight.seq <-data.frame(weight =seq(from =5, to =65, length.out =1e2) )# on récupère les prédictions du modèle pour ces valeurs de poidsmu <-data.frame(fitted(mod8, newdata = weight.seq, probs =c(0.055, 0.945) ) ) %>%bind_cols(weight.seq)pred_height <-data.frame(predict(mod8, newdata = weight.seq, probs =c(0.055, 0.945) ) ) %>%bind_cols(weight.seq)# on affiche les 6 premières lignes de pred_heighthead(pred_height)
Bürkner, P.-C. (2017). brms : An R Package for Bayesian Multilevel Models Using Stan. Journal of Statistical Software, 80(1). https://doi.org/10.18637/jss.v080.i01
Gelman, A., & Pardoe, I. (2006). Bayesian measures of explained variance and pooling in multilevel (hierarchical) models. Technometrics, 48(2), 241–251. https://doi.org/10.1198/004017005000000517
Marsman, M., & Wagenmakers, E.-J. (2017). Three Insights from a Bayesian Interpretation of the One-Sided P Value. Educational and Psychological Measurement, 77(3), 529–539. https://doi.org/10.1177/0013164416669201
Marsman, M., Waldorp, L., Dablander, F., & Wagenmakers, E.-J. (2019). Bayesian estimation of explained variance in ANOVA designs. Statistica Neerlandica, 0(0), 1–22. https://doi.org/10.1111/stan.12173
